Analysis of a typical HDX-MS experiment
We have seen the basic aspects of our functional modelling approach. We now
wish to roll out our method across all peptides in the experiment. The
fitUptakeKinetics function allows us to apply our modelling approach across
all the peptide in the experiment. We need to provide a QFeatures object
and the features for which we are fitting the model. The design will be extracted
from the column names or you can provide a design yourself. The parameter
initilisation should also be provided. Sometimes the model can’t be fit on the
kinetics. This is either because there is not enough data or through lack of
convergence. An error will be reported in these cases but this should not
perturb the user. You may wish to try a few starting values if there
excessive models that fail fitting.
res <- fitUptakeKinetics(object = MBPqDF[,c(1:24)],
feature = rownames(MBPqDF[,c(1:24)])[[1]],
start = list(a = NULL, b = 0.001, d = NULL, p = 1))
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
The code chunk above returns a class HdxStatModels indicating that a number
of models for peptide have been fit. This is simply a holder for a list
of HdxStatModel instances.
res
## Object of class "HdxStatModels"
## Number of models 102
We can easily examine indivual fits by going to the underyling HdxStatModel
class:
res@statmodels[[1]]@vis + scale_color_manual(values = brewer.pal(n = 2, name = "Set2"))
## Warning in brewer.pal(n = 2, name = "Set2"): minimal value for n is 3, returning requested palette with 3 different levels
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
We now wish to apply statistical analysis to these fitted curves. Our approach
is an empirical Bayes testing procedure, which borrows information across peptides
to stablise variance estimates. Here, we need to provide the original data
that was analysed and the HdxStatModels class. The following code chunk
returns an object of class HdxStatRes. This object tell us that statistical
analysis was performed using our Functional model.
out <- processFunctional(object = MBPqDF[,1:24], params = res)
## Warning in stats::optim(x = c(0.00109178357660077, 4.02437669894543e-05, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in densfun(x, parm, ...): NaNs produced
## Warning in densfun(x, parm, ...): NaNs produced
## Warning in densfun(x, parm, ...): NaNs produced
## Warning in stats::optim(x = c(0.00153823834115607, 0.00188329333146131, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in densfun(x, parm, ...): NaNs produced
out
## Object of class "HdxStatRes"
## Analysed using Functional model
The main slot of interest is the results slot which returns quantities of
interest such as p-values and fdr corrected p-values because of multiple testing.
The following is the DataFrame of interest.
out@results
## DataFrame with 102 rows and 9 columns
## Fstat.Fstat Fstat.numerator Fstat.denomenator
## <list> <list> <list>
## VIWINGDKGYNG_2 0.709762 0.00109178 0.00153824
## VIWINGDKGYNGL_2 0.0213688 4.02438e-05 0.00188329
## GDKGYNGLAEVG_3 0.980103 0.00129713 0.00132346
## LAEVGKKFEKDTGIKVTVEHPDK_4 0.809885 0.00117025 0.00144496
## AEVGKKFEKDTGIKV_4 0.806449 0.00140588 0.00174329
## ... ... ... ...
## YAVRTAVINAASGRQTVDE_3 -3.96761 -2.79196 0.703687
## AVINA_1 -3.87726 -0.213395 0.0550377
## AVINAASGRQTVDEAL_2 0.905874 0.00277575 0.00306417
## ASGRQTVDE_2 0.37392 0.00288595 0.00771809
## ASGRQTVDEAL_2 0.364883 0.00215214 0.00589818
## pvals fdr ebayes.pvals ebayes.fdr
## <numeric> <numeric> <numeric> <numeric>
## VIWINGDKGYNG_2 0.597056 1 0.589127 1
## VIWINGDKGYNGL_2 0.999008 1 0.998946 1
## GDKGYNGLAEVG_3 0.445918 1 0.442886 1
## LAEVGKKFEKDTGIKVTVEHPDK_4 0.536935 1 0.530748 1
## AEVGKKFEKDTGIKV_4 0.538920 1 0.526604 1
## ... ... ... ... ...
## YAVRTAVINAASGRQTVDE_3 1.000000 1 1.000000 1
## AVINA_1 1.000000 1 1.000000 1
## AVINAASGRQTVDEAL_2 0.483850 1 0.457787 1
## ASGRQTVDE_2 0.823819 1 0.803136 1
## ASGRQTVDEAL_2 0.830029 1 0.811177 1
## empirical.fdr fitcomplete
## <numeric> <integer>
## VIWINGDKGYNG_2 0.727165 1
## VIWINGDKGYNGL_2 0.762892 2
## GDKGYNGLAEVG_3 0.727165 3
## LAEVGKKFEKDTGIKVTVEHPDK_4 0.727165 4
## AEVGKKFEKDTGIKV_4 0.727165 5
## ... ... ...
## YAVRTAVINAASGRQTVDE_3 1.000000 98
## AVINA_1 1.000000 99
## AVINAASGRQTVDEAL_2 0.727165 100
## ASGRQTVDE_2 0.727165 101
## ASGRQTVDEAL_2 0.727165 102
We can now examine the peptides for which the false discovery rate is less
than 0.05
which(out@results$ebayes.fdr < 0.05)
## EAAFNKGETAM_2 AAFNKGETA_2 AAFNKGETAM_2 WYAVRTAVINA_2
## 55 56 57 96
Let us visualise some of these examples:
res@statmodels[[42]]@vis + res@statmodels[[45]]@vis
As we can see our model has picked up some subtle differences, we can further
visualise these using a forest plot. We can see the the functions are very similar
as the parameters are almost identical (a,b,p,d). However, we can see that
the deuterium differences are lower in 10% structural variant condition.
fp <- forestPlot(params = res@statmodels[[42]])
We can produce a table to actual numbers. We see that at all 4 timepoints
the deuterium difference is negative, though the confidence intervals overlap
with 0. Our functional approach is picking up this small but reproducible difference.
knitr::kable(fp$data)
It is also possible to visualize, these plots on a different scale. Of course,
changing the natural scaling will emphasis different parts of the plot and
could distort interpretation. In particular, if a log transform is used then
care should be taken when interpreting values around 0. We suggest examining
the numerical values in a forest plot or table alongside any transformation of
the variables. We suggest using the pseudo log transform as this allows
control the linearity of the plot, clearly demonstrating this a choice
of visualisation (and not of statistical modelling). The parameter sigma
below controls the scaling factor of the linear part of the transformation.
res@statmodels[[42]]@vis + scale_x_continuous(
trans = pseudo_log_trans(base = 10, sigma = 0.01), breaks = c(0, 10^(1:7)))

res@statmodels[[42]]@vis + scale_x_continuous(
trans = pseudo_log_trans(base = 10, sigma = 0.0001), breaks = c(0, 10^(1:7)))

res@statmodels[[42]]@vis + scale_x_continuous(
trans = pseudo_log_trans(base = 10, sigma = 10), breaks = c(0, 10^(1:7)))

Let’s us now have a look a situation where the changes are more dramatic.
res_wt <- fitUptakeKinetics(object = MBPqDF[, c(61:100)],
feature = rownames(MBPqDF[, c(61:100)])[[1]],
start = list(a = NULL, b = 0.001, d = NULL, p = 1))
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## Error in nlsModel(formula, mf, start, wts) :
## singular gradient matrix at initial parameter estimates
## [1] "Could not fit model, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
out_wt <- processFunctional(object = MBPqDF[, c(61:100)], params = res_wt)
## Warning in stats::optim(x = c(0.0840870774115754, 0.0705385629222093, 0.0354279988095238, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in stats::optim(x = c(0.00199791426867962, 0.00247716448153206, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in densfun(x, parm, ...): NaNs produced
We can visualise some of the result and generate plots.
res_wt@statmodels[[27]]@vis/res_wt@statmodels[[28]]@vis + plot_layout(guides = "collect")|(forestPlot(params = res_wt@statmodels[[27]], condition = c("WT", "W169G"))/forestPlot(params = res_wt@statmodels[[28]], condition = c("WT", "W169G")) + plot_layout(guides = "collect")) +
plot_annotation(tag_levels = 'a') + plot_layout(widths = c(1, 1))

An epitope mapping experiment
We now describe the analysis of an epitope mapping experiment. Here, the data
analysis is more challenging, since only 1 replicate in each condition, apo and
antibody, was performed. If we make some simplifying assumptions rigorous
statistical analysis can still be performed.
The experiment was performed on HOIP-RBR, we loaded the data below from inside
the package
HOIPpath <- system.file("extdata", "N64184_1a2_state.csv", package = "hdxstats")
HOIP <- read.csv(HOIPpath)
unique(HOIP$State)
## [1] "apo" "dAb13_1" "dAb13_2" "dAb25_1" "dAb25_2" "dAb27_1" "dAb27_2"
## [8] "dAb2_1" "dAb2_2" "dAb6_1" "dAb6_2"
HOIP$Exposure <- HOIP$Exposure * 60 #convert to seconds
filter(HOIP, Sequence == unique(HOIP$Sequence[1])) %>%
ggplot(aes(x = Exposure,
y = Center,
color = factor(State, unique(HOIP$State)))) +
theme_classic() + geom_point(size = 3) +
scale_color_manual(values = colorRampPalette(brewer.pal(8, name = "Set2"))(11)) +
labs(color = "experiment", x = "Deuterium Exposure", y = "Deuterium incoperation")

As before we need to convert data to an object of classes QFeatures
for ease of analysis.
First, we put the data into a DataFrame object. Currently, its in long format
so we switch to a wide format
HOIP_wide <- pivot_wider(data.frame(HOIP),
values_from = Center,
names_from = c("Exposure", "State"),
id_cols = c("Sequence"))
Now remove all columns with only NAs
HOIP_wide <- HOIP_wide[, colSums(is.na(HOIP_wide)) != nrow(HOIP_wide)]
The colanmes are not very informative, provide in the format X(time)rep(repliate)cond(condition)
colnames(HOIP_wide)[-c(1)]
## [1] "0_apo" "30_apo" "300_apo" "0_dAb13_1" "30_dAb13_1"
## [6] "300_dAb13_1" "0_dAb13_2" "30_dAb13_2" "300_dAb13_2" "0_dAb25_1"
## [11] "30_dAb25_1" "300_dAb25_1" "0_dAb25_2" "30_dAb25_2" "300_dAb25_2"
## [16] "0_dAb27_1" "30_dAb27_1" "300_dAb27_1" "0_dAb27_2" "30_dAb27_2"
## [21] "300_dAb27_2" "0_dAb2_1" "30_dAb2_1" "300_dAb2_1" "0_dAb2_2"
## [26] "30_dAb2_2" "300_dAb2_2" "0_dAb6_1" "30_dAb6_1" "300_dAb6_1"
## [31] "0_dAb6_2" "30_dAb6_2" "300_dAb6_2"
new.colnames <- gsub("0_", "0rep1", paste0("X", colnames(HOIP_wide)[-c(1)]))
new.colnames <- gsub("rep1", "rep1cond", new.colnames)
# remove annoying % signs
new.colnames <- gsub("%", "", new.colnames)
# remove space (NULL could get confusing later and WT is clear)
new.colnames <- gsub(" .*", "", new.colnames)
Now, we can provide rownames and convert the data to a QFeatures object:
qDF <- parseDeutData(object = DataFrame(HOIP_wide),
design = new.colnames,
quantcol = 2:34,
rownames = HOIP_wide$Sequence)
As before, we can produce a heatmap, we perform a simple normalisation for
ease of visualisation:
mat <- assay(qDF)
mat <- apply(mat, 2, function(x) x - assay(qDF)[,1])
pheatmap(t(mat),
cluster_rows = FALSE,
cluster_cols = FALSE,
color = brewer.pal(n = 9, name = "BuPu"),
main = "HOIP RBR heatmap",
fontsize = 14,
legend_breaks = c(0, 2, 4, 6,8,10,12, max(assay(qDF))),
legend_labels = c("0", "2", "4", "6", "8","10", "12", "Incorporation"))

Let us first perform a quick test:
res <- differentialUptakeKinetics(object = qDF[,1:33],
feature = rownames(qDF)[[1]][3],
start = list(a = NULL, b = 0.01, d = NULL),
formula = value ~ a * (1 - exp(-b*(timepoint))) + d)
## Warning in differentialUptakeKinetics(object = qDF[, 1:33], feature =
## rownames(qDF)[[1]][3], : NAs introduced by coercion
res@vis+ scale_color_manual(values = colorRampPalette(brewer.pal(8, name = "Set2"))(11))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

Whilst this analysis performs good fits for the functions, there are too many
degrees of freedom to perform sound statistical analysis. Hence, we normalize
to remove the degree of freedom for the intercept. For simplicity and to preserve
the original matrix, we reprocess the data. We then fit a simplified kinetic
model, where only the plateau is inferred.
cn <- new.colnames[c(1:3,10:12)]
HOIP_wide_nrm <- data.frame(HOIP_wide)
HOIP_wide_nrm[, c(2:4)] <- HOIP_wide_nrm[,c(2:4)] - HOIP_wide_nrm[,c(2)] # normalise by intercept
HOIP_wide_nrm[, c(11:13)] <- HOIP_wide_nrm[,c(11:13)] - HOIP_wide_nrm[,c(11)] # normalised by intercept
newqDF <- parseDeutData(object = DataFrame(HOIP_wide_nrm),
design = cn,
quantcol = c(2:4, 11:13), rownames = HOIP_wide$Sequence)
res_all <- fitUptakeKinetics(object = newqDF[,1:6],
feature = rownames(newqDF[,1:6])[[1]],
start = list(a = NULL),
formula = value ~ a * (1 - exp(-0.05*(timepoint))), maxAttempts = 1)
funresdAb25_1 <- processFunctional(object = newqDF[,1:6],
params = res_all)
## Warning in stats::optim(x = c(0.0236132769138663, 0.0111598488727292, 0.130743883316999, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in stats::optim(x = c(0.0466717805074177, 0.0136102073076068, 0.041005949389079, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in densfun(x, parm, ...): NaNs produced
We can have a look at the results:
funresdAb25_1@results
## DataFrame with 110 rows and 9 columns
## Fstat.Fstat Fstat.numerator Fstat.denomenator
## <list> <list> <list>
## GPGQECA 0.505943 0.0236133 0.0466718
## CAVCGWALPHNRM 0.819962 0.0111598 0.0136102
## CAVCGWALPHNRMQAL 3.18841 0.130744 0.0410059
## CAVCGWALPHNRMQALTSCE 1.51203 0.598316 0.395703
## AVCGWALPHNRM 0.138914 0.00401256 0.0288852
## ... ... ... ...
## ATERYLHVRPQPLAGEDPPAYQARL 0.882916 0.0634907 0.0719102
## ERYLHVRPQPLAGEDPPAYQ 0.841277 0.0671491 0.0798181
## ERYLHVRPQPLAGEDPPAYQARL 2.11026 0.172771 0.081872
## RYLHVRPQPLAGEDPPAYQARL 1.65969 0.120275 0.0724686
## LQKLTEEVPLGQSIPRRRK 4.48223 0.196504 0.0438406
## pvals fdr ebayes.pvals ebayes.fdr
## <numeric> <numeric> <numeric> <numeric>
## GPGQECA 0.516181 0.652643 0.479007 0.612683
## CAVCGWALPHNRM 0.416403 0.558589 0.406959 0.545920
## CAVCGWALPHNRMQAL 0.148709 0.320744 0.123017 0.300707
## CAVCGWALPHNRMQALTSCE 0.286210 0.499732 0.237829 0.441045
## AVCGWALPHNRM 0.728272 0.801099 0.708845 0.779730
## ... ... ... ... ...
## ATERYLHVRPQPLAGEDPPAYQARL 0.400605 0.557804 0.3563859 0.502596
## ERYLHVRPQPLAGEDPPAYQ 0.410930 0.558053 0.3658748 0.509446
## ERYLHVRPQPLAGEDPPAYQARL 0.219967 0.443454 0.1817334 0.386882
## RYLHVRPQPLAGEDPPAYQARL 0.267115 0.481684 0.2263725 0.436859
## LQKLTEEVPLGQSIPRRRK 0.101671 0.266280 0.0814249 0.246599
## empirical.fdr fitcomplete
## <numeric> <integer>
## GPGQECA 0.795215 1
## CAVCGWALPHNRM 0.778284 2
## CAVCGWALPHNRMQAL 0.778284 3
## CAVCGWALPHNRMQALTSCE 0.778284 4
## AVCGWALPHNRM 0.816964 5
## ... ... ...
## ATERYLHVRPQPLAGEDPPAYQARL 0.778284 106
## ERYLHVRPQPLAGEDPPAYQ 0.778284 107
## ERYLHVRPQPLAGEDPPAYQARL 0.778284 108
## RYLHVRPQPLAGEDPPAYQARL 0.778284 109
## LQKLTEEVPLGQSIPRRRK 0.778284 110
which(funresdAb25_1@results$ebayes.fdr < 0.05)
## IQLRESLEPDA RESLEPDAYALFHKKLTEGVL YALFHKKLTEGVL
## 36 42 43
## REQLEATCPQCHQTF EATCPQCHQTF MYLQENGIDCPKCKF
## 52 53 65
## YLQENGIDCPKCKFSYA LQENGIDCPKCKFSYA
## 68 70
We can plot these kinetics to see what is happening. This allows us to visualise
region of protection and deprotection, potentially identifiying the epitope.
(res_all@statmodels[[36]]@vis +
res_all@statmodels[[42]]@vis +
res_all@statmodels[[43]]@vis +
res_all@statmodels[[65]]@vis +
res_all@statmodels[[68]]@vis +
res_all@statmodels[[70]]@vis +
res_all@statmodels[[52]]@vis +
res_all@statmodels[[53]]@vis ) + plot_layout(guides = 'collect')
We can make a Manhattan plot to better specially visualise what’s happening.
#We need to provide an indication of "difference" so we can examine deprotected
# or prected regions
diffdata <- assay(newqDF)[,6] - assay(newqDF)[,3]
sigplots <- manhattanplot(params = funresdAb25_1,
sequences = HOIP$Sequence,
region = HOIP[, c("Start", "End")],
difference = diffdata,
nrow = 1)
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
sigplots[[1]] + plot_layout(guides = 'collect')

We can visualise this in a peptide plot which helps us understand the nature
of the overlap
fpath <- system.file("extdata", "HOIP.txt", package = "hdxstats", mustWork = TRUE)
HOIPfasta <- readAAStringSet(filepath = fpath, "fasta")
scores <- funresdAb25_1@results$ebayes.fdr
out <- plotEpitopeMap(AAString = HOIPfasta[[1]],
peptideSeqs = unique(HOIP$Sequence),
numlines = 2,
maxmismatch = 1,
by = 1,
scores = 1 * (-log10(scores[unique(HOIP$Sequence)]) > -log10(0.05)) + 0.0001,
name = "significant")
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
out[[1]]/(out[[2]]) + plot_layout(guides = 'collect') & theme(legend.position = "right")

We can further visualise this a barcode of particular residues, here we use
residue level averaging to obtain results at the residue level.
scores <- funresdAb25_1@results$ebayes.fdr
out2 <- plotEpitopeMapResidue(AAString = HOIPfasta[[1]],
peptideSeqs = unique(HOIP$Sequence),
numlines = 2,
maxmismatch = 1,
by = 5,
scores = scores[unique(HOIP$Sequence)],
name = "-log10 p value")
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
out2[[1]]/out2[[2]] + plot_layout(guides = 'collect') & theme(legend.position = "right")
We can also plot multiple residue maps on the same plot so that we can compare
different antibodies.
scores <- funresdAb25_1@results$ebayes.fdr
avMap25_1 <- ComputeAverageMap(AAString = HOIPfasta[[1]],
peptideSeqs = unique(HOIP$Sequence),
numlines = 2, maxmismatch = 1,
by = 10, scores = scores[unique(HOIP$Sequence)],
name = "-log10 p value")
## generate results from other dAB
cn <- new.colnames[c(1:3,19:21)]
HOIP_wide_nrm <- data.frame(HOIP_wide)
HOIP_wide_nrm[,c(2:4)] <- HOIP_wide_nrm[,c(2:4)] - HOIP_wide_nrm[,c(2)]
HOIP_wide_nrm[,c(20:22)] <- HOIP_wide_nrm[,c(20:22)] - HOIP_wide_nrm[,c(20)]
newqDF2 <- parseDeutData(object = DataFrame(HOIP_wide_nrm),
design = cn,
quantcol = c(2:4,20:22),
rownames = HOIP_wide$Sequence)
res_all2 <- fitUptakeKinetics(object = newqDF2[,1:6],
feature = rownames(newqDF2[,1:6])[[1]],
start = list(a = NULL),
formula = value ~ a * (1 - exp(-0.07*(timepoint))), maxAttempts = 1)
## Warning in max(data$value): no non-missing arguments to max; returning -Inf
## Warning in differentialUptakeKinetics(object = object, feature = x, start =
## start, : NAs introduced by coercion
## [1] "model fit failed, likely exessive missing values"
funresdAb27_2 <- processFunctional(object = newqDF[,1:6],
params = res_all2)
## Warning in stats::optim(x = c(9.29818825569032e-05, 0.090283790910004, 0.166406795432091, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in densfun(x, parm, ...): NaNs produced
## Warning in stats::optim(x = c(0.00297145686513731, 0.032686171739902, 0.0181780887330962, : one-dimensional optimization by Nelder-Mead is unreliable:
## use "Brent" or optimize() directly
## Warning in densfun(x, parm, ...): NaNs produced
scores <- funresdAb27_2@results$ebayes.fdr
# compute average map
avMap27_2 <- ComputeAverageMap(AAString = HOIPfasta[[1]],
peptideSeqs = unique(HOIP$Sequence),
numlines = 2,
maxmismatch = 1,
by = 10,
scores = scores[unique(HOIP$Sequence)],
name = "-log10 p value")
# set rownames
rownames(avMap25_1) <- "dAb25_1"
rownames(avMap27_2) <- "dAb27_2"
# store in a list
avMap <- list(avMap27_2 = avMap27_2,
avMap25_1 = avMap25_1)
#plotting
out3 <- plotAverageMaps(avMap, by = 20)
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
out3[[1]]/out3[[2]] + plot_layout(guides = 'collect') & theme(legend.position = "right")
